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We introduce a simple spherical model whose structural properties are similar to the ones gen- 
erated by models with directional interactions, by employing a binary mixture of large and small 
hard spheres, with a square-well attraction acting only between particles of different size. The small 
particles provide the bonds between the large ones. With a proper choice of the interaction param- 
eters, as well as of the relative concentration of the two species, it is possible to control the effective 
valence. Here we focus on a specific choice of the parameters which favors tetrahedral ordering and 
study the equilibrium static properties of the system in a large window of densities and tempera- 
tures. Upon lowering the temperature we observe a progressive increase in local order, accompanied 
by the formation of a four-coordinated network of bonds. Three different density regions are ob- 
served: at low density the system phase separates into a gas and a liquid phase; at intermediate 
densities a network of fully bonded particles develops; at high densities — due to the competition 
between excluded volume and attractive interactions — the system forms a defective network. The 
very same behavior has been previously observed in numerical studies of non-spherical models for 
molecular liquids, such as water, and in models of patchy colloidal particles. Differently from these 
models, theoretical treatments devised for spherical potentials, e.g. integral equations and ideal 
mode coupling theory for the glass transition can be applied in the present case, opening the way 
for a deeper understanding of the thermodynamic and dynamic behavior of low valence molecules 
and particles. 

PACS numbers: 82.70.Gg, 82.70.Dd, 61.20.Lc 



I. INTRODUCTION 



The study of dynamic arrest in atomic and molecu- 
lar systems is an active field of research^ . Close to the 
glass transition, in a small window of values of the ex- 
ternal parameters (temperature or density/pressure) the 
dynamics of the system slows down by fifteen or more 
orders of magnitude. The slowing down of the dynamics 
varies in the temperature or density/pressure dependence 
from an Arrhenius law for the so-called strong network 
glass-forming liquids to a more complex super- Arrhenius 
behavior for fragile glass-forming liquids [2J. Signatures of 
the slowing down of the dynamics are observed at tem- 
peratures higher (or density lower) than the calorimetric 
glass transition and have been interpreted as the genuine 
precursor of the arrest phenomenon by the ideal Mode 
Coupling Theory (MCT) [3] . The study of the dynamics 
in colloidal systems [U [5] has added new fuel to the dis- 
cussion on dynamic arrest. Experimental realizations of 
simple models amenable to theoretical investigation has 
considerably boosted the understanding of the glass phe- 
nomenon and opened up an ampler view by adding to the 
arena the gelation issue, i.e. the possibility of observing 
arrest at low density, driven by the formation of inter- 
particle attractive bonds. Indeed, nowadays it is possi- 
ble to realize colloidal systems which closely follow the 
hard-sphere (HS) equation of state, as well as to control 
an additional attractive interaction, both in range and 
in strength, moving the field of colloidal liquids towards 
molecular systems [B]. 



Recent studies have shown that when the hard-core is 
complemented by a spherical attractive potential, phase 
separation preempts the possibility of continuously ap- 
proaching the slowing down of the dynamics at low den- 
sity. For Lennard-Jones particles, Sastry[7] showed that 
the glass line intersects the liquid-gas spinodal on the liq- 
uid branch, suggesting that homogeneous arrested states 
are only possible for significantly large densities. More re- 
cently, the same scenario has been shown to hold even in 
the limit of very short-ranged spherical attractions [HI [5], 
down to the Baxter limit (infinitesimal attraction range) . 
On the other hand, a series of studies have suggested 
that a progressive arrest at small packing fraction can 
be observed if the inter-particle interactions are highly 
directional (patchy interactions), when the effective va- 
lence becomes small 1 0J . On decreasing the valence below 
six, the gas-liquid unstable region progressively shrinks 
to smaller densities, opening up an intermediate region 
where a stable network of bonded particles forms. The 
shrinking of the unstable region can be tuned continu- 
ously down to a vanishing densities on decreasing the 
valence down to two [IT]. When the average valence is 
slightly larger than two, it is possible to observe empty 
liquids, i.e. states with a temperature smaller than the 
critical temperature but with an extremely small liquid- 
state density. Numerical studies of the dynamics of liq- 
uids with patchy interactions have shown that the bond 
lifetime of the liquid network progressively increases upon 
cooling, providing evidence of the possibility of approach- 
ing arrested states continuously. The slowing down of the 
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dynamics in the newly accessible density region follows an 
Arrhenius Iaw|12 [ fl3 | ri4]. The arrest transition can be si- 
multaneously interpreted in terms of equilibrium gelation 
in the colloid community [51 H2"] and of strong- glass forma- 
tion in the supercooled liquids community 1 1 3 1 [T^] . These 
studies have clearly shown that reduction of valence is an 
essential ingredient for extending the glass line from the 
high-density/pressure region, dominated by repulsion, to 
intermediate regions where attraction and repulsion co- 
operate, down to the region of lower densities where at- 
traction become the dominant arrest mechanism. 

While simulations are a well documented method for 
studying the static and dynamic properties of non- 
spherical models, theoretical approaches are still not well 
developed, especially when the geometry of the bonds is 
such that inter-particle correlations propagate over sev- 
eral neighbors. The situation is even worse for micro- 
scopic theories of the glass transition, requiring struc- 
tural quantities as input data. For the case of molecules, 
molecular [HI HU ED EB] and site-site PH H3 125 exten- 
sions of MCT have been developed, but their use has 
been limited due to the nature of the approximations 
and to the complexity of the calculations. An effective 
spherical approximation for non-spherical potentials has 
also been recently reported [25]. For these reasons, the 
theoretical evaluation of the MCT glass line for patchy 
non-spherical models has never been attempted so far. 
For the case of water (at ambient pressure) a molecular- 
MCT calculation has been reported |16l I23j . Besides wa- 
ter, the only network forming system whose slow dynamic 
properties have been investigated — starting from ap- 
propriate structural properties — is silica [2U [23] [25]. In 
the investigated model of silica, the network is formed 
by a binary mixture of two spherically interacting parti- 
cles and hence all the complications associated to angular 
constraints are missing. Theoretical studies of silica have 
been limited to one specific value of the density and no 
clear picture of the glass-line in the T-p plane has been 
provided. 

As a first step toward a simple network-forming model 
for which theoretical calculations of the location of the 
arrest line in the phase diagram are in principle feasible, 
we introduce here a simple binary mixture model of par- 
ticles interacting via short-range square well potentials 
and which is able to reproduce the same pattern (both 
for structural and dynamic quantities) of the previously 
studied non-spherical patchy models. The idea is bor- 
rowed from existing models for silica based on pair-wise 
additive interactions [27 but with a much simpler inter- 
action potential. Due to the absence of the long-range 
electrostatic interactions, the model can be simulated in 
the entire phase diagram and the slowing down of the 
dynamics can be followed over a window of more than 
five orders of magnitude. The present article provides 
an evaluation of the structural properties of this model, 
based on extensive event-driven molecular dynamics sim- 
ulations. A future companion article [28 will report the 
dynamic properties of the model. The manuscript is or- 



ganized as follows: section II introduces the model; sec- 
tion III reports results for the structural, energetic and 
geometric properties . Section IV discusses the phase 
diagram, complemented by the study of various thermo- 
dynamic loci, while section V is devoted to conclusions. 



II. THE MODEL 

Previously studied simple models for network forming 
liquids are based on many-body interactions [TU1 [2§] or 
angular constraints [TT] [30]. Here we introduce a simple 
model retaining both pair-wise additivity and spherical 
interactions, but which is capable of producing geomet- 
rical arrangement into a locally ordered network struc- 
ture. The oxymoron spherical model with directional in- 
teractions is realised through the introduction of a second 
species in the mixture that represents a sort of floating 
bond. In this way, an effective one-component directional 
potential for the colloidal particles is obtained. 

In more details, we consider N\ hard sphere colloids 
(for which we reserve the name ■particles in the follow- 
ing) with diameter an. We mix them with N 2 small 
particles with hard-sphere diameter 1722 (named floating 
bonds in the following). Particles and floating bonds in- 
teract via a non-additive (NA), short-ranged, attractive 
square-well (SW) potential of depth — uq and range S. 
Thus, the floating bonds link particles providing a con- 
nection between them. More precisely the interaction 
potential is 
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We choose the potential parameters in such a way that 
(i) each floating bond binds no more than two particles; 
and (ii) the maximum number of floating bonds bind- 
ing to a single particle is fixed, providing the valence of 
the model. To enforce condition (i), we recall that three 
identical touching hard-spheres of diameter an create a 
cavity which can incorporate a hard-sphere with diam- 
eter up to d c = (2/-\/(3) — l)<7n = 0.1547(Tii making 
contact with each of the three spheres. Hence, if the ge- 
ometric condition a^ 2 A + 5 < (an + d c )/2 is met, the 
floating bond can only be simultaneously involved in two 
attractive interactions. We choose a^ 2 A = 0.55cth and 
5 = 0.03cr^ A so that 2(a^ A + S) - an = 0.134cr n < d c . 
As a result of this choice, a floating bond can be isolated 
(with potential energy zero), bonded to one particle (po- 
tential energy — uq) or bonded to two distinct particles 
(potential energy —2u ). Only in the last case, the float- 
ing bond provides a link between the two particles, which 
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are thus considered bonded and belonging to the same 
cluster. 

To satisfy the limited valence condition (ii) we fix the 
hard-sphere diameter of the floating bonds (722, deter- 
mining the closest distance between them. To model 
a system with valence four (with tetrahedral ordering) 
we choose 022 = 0.8<7n (and hence, the mixture con- 
sidered here has a negative non-additive parameter A = 
cr i2 j4 /( cr ii + °22) — 1 — —0.69). This choice is dictated by 
geometric considerations. Indeed, the distance between 
vertices of a perfect tetrahedron is 2V6/3 times the dis- 
tance between the center and the vertex. The distance of 
an interacting floating bond from the center of a particle 
varies between af 2 A and a^ 2 A + S. Hence, a tetrahedral 
arrangement requires a value of 0-22 < 0.93<7n. To pre- 
vent the possibility of a planar square arrangement of 
bonds (a geometry which would allow a valence of six, 
with four bonds on the equatorial plane and two bonds 
on the poles) we need to impose 022 > 0.802<th. We 
have checked that our choice 022 = 0.8<7n never gener- 
ates particles connected to more than four floating bonds. 
Finally we fix the number ratio between particles and 
floating bonds following the stoichiometry of the system, 
i.e. by imposing that in the fully bonded ground state all 
particles participate to four bonds and that all floating 
bonds have energy —2u ., i.e. N 2 /Ni = 2. 

We study a system of Ni = 1000 particles and per- 
form standard event-driven MD simulations in the NVE 
ensemble. For the smallest studied temperatures, equi- 
libration requires several months of CPU time. Packing 
fraction is defined as (j> = irNi(an/ L) 3 /6, i.e. as the 
fraction of volume occupied only by the particles. In 
these units, taking into account the minimal and max- 
imal bond distance between two particles, the diamond 
crystal is mechanically stable between 0.233 < <j> < 0.255 
(since a diamond crystal of touching hard-spheres has 
4> — 0.340). Units of length and energy are an and uq, 
while the Boltzmann constant ks is set equal to 1. We 
equilibrated the mixture for a wide range of <p an d T, 
up to = 0.56 and down to T = 0.065. We notice that 
only for <p > 0.52 crystallization (into a fee structure) 
was sometimes (but not always) detected at low T. 

The use of a square-well potential to model the inter- 
actions makes it possible to unambiguously define the 
existence of a bond between two large particle. Indeed, 
if the energy of a floating bonding particle is —2uq, then 
the two large particles interacting with the floating bond 
particle are considered bonded. In this way, large par- 
ticles can be partitioned into different clusters and the 
connectivity properties of them can be examined. More 
explicitly, to test for percolation, the simulation box is 
duplicated in all directions, and the ability of the largest 
cluster to span the replicated system is controlled. If the 
cluster in the simulation box does not connect with its 
copy in the duplicated system, then the configuration is 
assumed to be nonpercolating. The boundary between a 
percolating and a nonpercolating state point is then de- 
fined as the probability of observing infinite clusters in 




FIG. 1: Partial radial distribution function of particles 311(7") 
along a \ow-(f> isochore, i.e. <j) = 0.25, and upon decreasing T. 

50% of the configurations. 



III. STATIC PROPERTIES 
A. Radial distribution functions 

To visualize how directional interactions build up for 
large particles upon decreasing temperature we report 
the behaviour of the particles radial distribution function 
<7ii(r) in Fig. [I] At large T, no significant correlations 
are present and the system behaves almost as a hard- 
sphere fluid mixture. The first peak of <?ii(r) is found at 
(Tix. However, as T is lowered, more and more particles 
are bonded within the attractive well, so that for <7n < 
r < 2(af i 2 A + 8) a larger and larger correlation is built up. 
For r > 2{af 2 A + S) an anti-correlation develops. When 
T < 0.1, the peak of gn( r ) shifts from cr n to 2a^ 2 A , a 
signature of the progressive role played by the floating 
bonds in structuring the particles: while the growth at 
an saturates, that at 2af 2 A continues to increase. 

A clear signal of the progressive structuring upon cool- 
ing is provided by the increase of the nearest neigh- 
bors peaks. The second peak is centered at a distance 
~ 1.78(7ii, while the third one is found at ~ 2.53ern. 
Such a sequence of peak positions, with non-integer ra- 
tios, is typical of network-forming systems with tetra- 
hedral arrangement [3TJ [351 [33]. It is also important to 
notice that, for T < 0.09, such peaks do not show a signif- 
icant variation in intensity, indicating that the structur- 
ing process is essentially completed and that the system 
has approached an almost perfect long-range tetrahedral 
arrangement. 

The density dependence of the partial radial distribu- 
tion functions is reported, for T — 0.09, in Fig. [2j We 
first focus on the evolution of (?n(r) in the main panel. 
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FIG. 2: Partial radial distribution functions of particles gn(r) 
at fixed low T — 0.09 and varying (f>. In the insets also partial 
radial distribution functions of floating bonds 322(f) and of 
mixed type (712(f) are shown. 



Data show that the tetrahedral order, clearly visible at 
small and intermediate densities (0.225 < <f> < 0.40), is 
progressively lost with increasing <f>. At this T, the curve 
corresponding to <p = 0.30 already shows small deviations 
in the location of the second peak. For larger <p the loca- 
tion of the second peak moves to smaller distances and its 
amplitude decreases. Now the first peak is always found 
at an rather than a^ 2 A for <j) > 0.30. A similar shift 
and rearrangement is observed for secondary peaks, as 
well as for the evolution of the other partial distribution 
functions 312 (f) and 322(f), shown in the insets. 

We further notice (not shown) that the 0-value where 
tetrahedricity starts to get progressively lost is found to 
depend on T. On further decreasing T, the maximum <fi 
with dominant local order increases: at T = 0.08, 0.075, 
also 4> = 0-35 becomes more tetrahedral- like. However, 
for tj> — 0.40, the structure even at the lowest equilibrated 
T, which is actually the closest point to the ground state 
structure that we could determine in our long-running 
simulations (see below), remains always far from that of 
a tetrahedral network. 

These results suggest that geometric frustration due 
to packing acts against the formation of a tetrahedral 
network above a certain <f>, giving rise to a competition 
between energetic (directional attraction) and entropic 
(excluded volume) interactions. 



B. Static structure factor 

We also study the behaviour of the nor- 
malized partial static structure factors, i.e. 
Sij(q) = (Ip^qJpiCqJD/WJV,-) 1 / 2 , where 
Pj'(q) = ^ m J =1 exp(iq'r m ) is the wave vector q 




qa 

FIG. 3: Normalized Sn(q) for <j> = 0.35 and varying T. Inset: 
Comparison with PY theoretical results. 



component of the density of species j (and the sum runs 
over all Nj particles of type j). We focus only on the 
partial structure factor of particles Sn(q), aiming to 
emphasize those features that are characteristic of the 
establishment of a fully connected network of bonds. 
Figure [3] shows the behavior of Sn(q) as a function of T 
at <j> = 0.35. We notice that, at high T, Sn(q) displays 
a nearest-neighbour peak centered around qa ~ 7, 
characteristic of the hard-sphere interactions and of all 
simple liquids. As T is decreased the height of such a 
peak also decreases, eventually giving rise at low T to 
a splitting in two peaks: one for larger length-scales at 
qa ~ 5 and one at smaller length-scales qa ~ 8.5. The 
first one seems to saturate at low T, while the second 
continues to increase in amplitude. The splitting of 
the main peak into two components is characteristic of 
network forming tetrahedral liquids |311 134) and it is 
associated to the formation of an open local structure. 

Simple integral equation theories are not able to pre- 
dict the splitting of the peak on cooling. We solved nu- 
merically the binary Ornstein-Zernike equation for our 
mixture within Percus-Yevick (PY) approximation [35 at 
the same packing fraction <p = 0.35 upon lowering T . Re- 
sults are shown in the inset of Fig. [3] together with a 
high-T numerical curve. While at high T PY provides 
a good description of the system, it fails severely upon 
lowering T, where the splitting of the main peak is not 
at all captured. Moreover, numerical convergence of the 
PY solution can not be achieved for T < 0.125. Hence, 
the PY solution is only able to capture a small decrease 
of the peak, without accounting for the formation of the 
tetrahedral network, as expected due to the symmetric 
approximation contained in PY. This stresses the im- 
portance to develop more elaborated integral equations 
which are able to account for the angular correlations in- 
troduced by the bonding. In this respect, a step forward 
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FIG. 4: Normalized Sii(g) for T = 0.09 and varying , 



could be made by exploring PY approximation applied 
to the associative Ornstein-Zernike equation (PY-AOZ) 
which generally provides better results, compared to sim- 
ulations, for network- forming liquids [36, 37 . 

The structure factor provides also information on the 
proximity to an unstable region since the low q behavior 
is related to the system compressibility. Indeed, close to 
an unstable state of the particles, Sn(q) increases signif- 
icantly at small q. Figure [4] shows 5*11(5) as a function 
of packing fraction, from (j> — 0.225 (just close to the 
phase separation liquid boundary) to <f> = 0.54 along a 
low temperature isotherm (T — 0.09). On increasing </>, 
the critical low q fluctuations disappear, giving rise to the 
network two-peak structure, then crossing continuously 
to a single-peaked standard Sn(q) that finally resembles 
that of high density simple liquids for <fi > 0.40. 

The study of the structural properties reported so far 
clearly indicates that the possibility of forming a well- 
connected disordered tetrahedral network arises only in 
a finite window of intermediate densities, i.e. 0.225 < 
4> It 0.40. For smaller </> values, particles are not close 
enough to form a fully bonded structure and phase sep- 
aration into a gas and a bonded liquid is preferred. For 
larger values, the local density around each particle be- 
comes more and more incompatible with the open struc- 
ture which is characteristic of tetrahedral networks and 
packing becomes the leading driving mechanism control- 
ling the structure. 



C. Energy 

In the present model, attraction takes place only be- 
tween particles and floating bonds. Moreover, the inter- 
action range is such that each floating bond can inter- 
act simultaneously only with up to two distinct parti- 
cles. This has the very important consequence that, for 



the chosen N1/N2 ratio, the energetic ground state of 
the system is known, being equal to two times the num- 
ber of floating bonds (in units of — Uq). In the present 
model, this means that a fully bonded configuration has 
a ground state energy E gs — —2uqN2 or equivalcntly 
E gs = ~^U[)Ni. We notice that an energetic gain is 
at hand both when a particle sticks to a floating bond 
(of contribution — u$) and when a true bond is estab- 
lished between two particles (of contribution — 2uq). To 
differentiate between these two situations, we use E to 
indicate the potential energy of the system (i.e. propor- 
tional the number of connections between a particle and 
a floating bond), while we use E^ to count the potential 
energy associated to particle-particle bond contributions 
(i.e. proportional to the number of floating bonds with 
potential energy — 2uo). Of course, the two quantities 
tend to become identical as T is lowered. We then define 
the bond probability as pi, = (Eb — E gs )/E gs . 

The possibility of knowing theoretically the ground 
state of the system (a property shared with other limited 
valence models [UJ US [TU ESI (33]) offers the possibil- 
ity of unambiguously tracking the approach to the fully 
bonded state and detect the range of densities where in- 
deed the system may reach continuously E gs . It has been 
previously suggested [3H] that the ability to reach in equi- 
librium almost fully bonded states is indeed a specific 
feature of network forming liquids. 

Fig. |5][a) shows the potential energy relative to the 
ground state energy, expressed as E/E gsi as a function 
of <p for all studied temperatures. A magnification in 
semi-logarithmic scale of the low-T isotherms is offered 
in Fig. [5jb) for the positive ratio (E — E gs )/E gs . Except 
for high temperatures, for which the <p dependence of E is 
monotonic (as in simple square- well models), a minimum 
in density appears. Since here the excluded volume inter- 
action is modeled via a hard-core interaction, the increase 
of E at large (f> can arise only from a progressive break- 
ing of the bonds. For intermediate T, when the average 
number of bonds per particle is still small, the minimum 
is met only at very large (f>. But, upon lowering T, bond- 
ing becomes more and more extensive and the density at 
the minimum progressively decreases, reaching <j) sa 0.3 
at the lowest studied T. The T-dependence of the density 
minimum is shown in the phase diagram reported below 
(see Fig. [10]). 

Fig. [5][b) also shows that, in a large window of den- 
sities, the system is able to essentially reach the disor- 
dered ground state. Indeed, more than 95% of the bonds 
(intended as E^) are formed. This window of densities 
is roughly the same for which structural properties sug- 
gested the presence of a tetrahedral network of bonds. 
Since particles have already formed almost all possible 
bonds at the lowest studied temperatures, a further de- 
crease in T, beyond the ones which we have been able to 
investigate, will not lead to further significant structural 
changes. 

Being the ground state energy a priori known, it is 
possible to investigate the T-dependence of the energy 
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FIG. 5: Density dependence of the potential energy per 
particle relative to the ground state energy along different 
isochores: (a) E/E ga for all studied state points and (b) 
(E — E gs ) /E gs in a magnification of the low T region in semi- 
logarithmic scale. 



on approaching the fully bonded state. To this aim, 
Fig. [6] shows E - E gs vs 1 JT in an Arrhenius plot for 
all investigated densities. Low <f> data are limited to the 
T-region above the liquid-gas instability. In the region 
where phase separation is not encountered, a clear Ar- 
rhenius dependence is found. However, for large enough 
densities, the energy decrease tends to saturate to a 
finite constant value. Only at intermediate densities 
(0.225 < 4> < 0.40), it appears to be possible to reach 
a fully bonded (defect free) network. At larger densities 
the lowest energy state appears to be larger than the fully 
bonded one, suggesting that the total volume constraint 
imposes the presence of a finite number of defects in the 
fully bonded network. This is not surprising in view of 
the fact that the establishment of a tetrahedrally coor- 
dinated environment enforces a constraint on the local 



FIG. 6: Potential energy per particle E — E gs vs. l/T. Note 
that at intermediate densities an Arrhenius approach to the 
ground state is observed. The dashed line is a reference curve 
with activation energy —0.5. 

density. 

To quantify the Arrhenius behaviour in an unbiased 
way, we fit the T-dependence of the potential energy with 
the functional form 

[E - E gs ] = E +B exp(E a /T) (2) 

in the region where a sufficient bonding is present. To re- 
alize this condition, we consider only state points where 
(E — E gs ) I Egs > 0.825, in order to have the same fitting 
conditions for all studied isochores. With the chosen fit- 
ting function, we can fully reproduce the behaviour of E 
both in the region where a fully connected network state 
is reached as well as that for states where frustration 
comes into play. Three fitting parameters are involved: 
.Bo, B and the activation energy E a . The results of the 
fits are summarized in Fig. [7] for all studied 4>. We ob- 
serve that the system approaches the expected ground 
state (Eq — 0) up to roughly ~ 0.40, while a consistent 
increase of E is found at larger tj>. Hence, beyond such 
4> value, the system can not reach a fully bonded config- 
uration due to the excluded-volume contributions in the 
free energy. Even at vanishing temperatures, when mini- 
mization of the energy is the dominant driving force, it is 
impossible for geometric reasons to reach a fully bonded 
state and the effective ground state, reached along a con- 
stant cj) cooling path, is characterized by a non-zero frac- 
tion of unformed bonds. Data in Fig. [6] show that it 
is possible to unambigously define an optimal region for 
network formation for 0.235 < <f> < 0.40 (the lower bound 
being fixed by the presence of phase separation) . In this 
density window, the fully bonded ground state can be 
reached in equilibrium on isochoric cooling. 

It is interesting to observe that the presence of a mini- 
mum in E(<p) along isotherms is consistent with a defect- 
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FIG. 7: Fit parameters from the fitting Arrhenius law in Eq. 
|2j effective ground-state energy ratio Eo/E ga (top) and acti- 
vation energy E a (bottom). 



free ground state, but only in a small region of densi- 
ties. The optimal network- forming region largely coin- 
cides with that where the structural properties also show 
tcthrahcdral arrangement, with the exception of the state 
points close to <fi = 0.40, where still E — > even though 
the network begins to be deformed due to the increase 
in packing. Such feature is only possible in the present 
model due to the flexibility of the chosen interaction pa- 
rameters. 

The ^-dependence of the activation energy (bottom 
panel) is also instructive. It is observed that E a is 
slightly larger than —0.5 within the optimal network- 
forming region, showing an almost monotonic decrease 
up to ~ 0.50, then a reversal of trend is observed at 
the highest studied values. Note that the value 0.5 is 
the theoretical value expected for the breaking process 
of independent bonds [30]. 

The low T Arrhenius behavior carries with it an- 
other important thermodynamic feature, the presence of 
a maximum in the constant volume V specific heat Cy. 
The behaviour of Cy along the studied isochores is re- 
ported in Fig. [8] In this respect, the present model con- 
firms the suggestion [39 that a line of Cy maxima in the 
phase diagram is present when bonding is the relevant 
driving force. Differently from the model in Ref. [39], 
a monotonic increase of the T- value of the Cy-maxima 
is observed upon increasing (f>. It is not a coincidence 
that maxima in Cy are also characteristic of reversible 
self-assembly processes |4 11 [42] . i.e. of systems in which 
bonding plays the leading role. The locus of Cy-maxima 
extracted from these results will be discussed in Sec. IIVI 
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FIG. 8: Constant volume specific heat Cy = (dE/dT)v for 
all studied isochores. 



D. Angular distribution of bonds 

A direct measure of the tetrahedral structure of the 
network is offered by the evaluation of the distribution 
of the angle 9 between bonded triplets of large particles, 
repeating the analysis that is usually performed for sim- 
ilar models j43] |44j [45] . The angle distribution P(9) is 
shown as a function of <f> for T = 0.075 in Fig. |9ja). 
For <j) < 0.30, P(0) is almost independent of <fi, suggest- 
ing the approach to a tetrahedral network with almost 
fully connected particles. A peak close to the tetrahe- 
dral angle is found, and an average value of the triplet 
angle of w 108.2. Also, around 8 = 0.60 a secondary 
peak is observed, a geometric signature of the presence 
of a small number of three-particle rings |43M45j . A Gaus- 
sian fit is attempted (dashed curve in Fig. ga) ), which 
works well only in the large-0 region, due to the presence 
of rings. In agreement with the structural indicators dis- 
cussed previously, P{6) deviates more and more from the 
typical tetrahedral shape on increasing <j>. When the ex- 
cluded volume contribution becomes dominant over the 
attractive interactions, the main tetrahedral peak is pro- 
gressively lost and a wider distribution is observed, ac- 
companied also by a smaller average angle and by an 
increase of three-particles loops. 

In Fig. [9jb) a comparison of P{9) for different 4- 
coordinated models is offered, completing the data al- 
ready presented in [14 . In that work, P{6) for a model 
without any geometric correlation in the location of the 
bonding sites (the N max model [12] [46] [47] ) — showing a 
rather flat distribution of angles due to the lack of pre- 
ferred orientation between the bonds (i.e. random bond 
organization) — was compared with P{9) for two mod- 
els with well defined locations of the bonding sites on the 
particles surface, the so-called primitive models of wa- 
ter (PMW)[48^ and primitive models of silica (PMS)[4"9"]. 
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FIG. 9: Distribution of the angle 9 between bonded particles 
triplets. Full thin lines are guides to the eye. (a) P(6) for 
various (f> at T = 0.075. The thick dashed line is a Gaussian 
fit to the low-0 states, suggesting an average angle ~ 107.4 
and variance ~ 29.0; (b) Comparison of P{9) in the optimal 
network-forming density [cj> — 0.25, T = 0.075- triangles ) 
with different 4-coordinated models: PWM (dot-dashed line), 
PSM (full line), 7V max (dashed line) (all from Ref. [13]). as well 
as for the present model at two different state points (triangles 
and circles). 



The PMW, at its optimal density, displays a quite sharp 
and narrow peak, with an average angle ~ 109 degrees, 
due to the rigid organization of the bonding sites model- 
ing the hydrogen atoms and the lone-pairs of the water 
molecule, with a roughly Gaussian shape of P(6). On the 
other hand, the PMS at its optimal network-forming den- 
sity has a wider and more asymmetric distribution, which 
is peaked at lower angles 6 « 90 degrees. In addition, the 
N m ax model was found to exhibit a dominant contribu- 
tion of three-particles loops or touching triplets (6 = 60 
degrees) , a feature that was absent in PMW and negligi- 
ble in PMS. Comparing the results for the present binary 



model to the existing data, we observe that, within the 
optimal network-forming region, the organization of par- 
ticles is very reminiscent of that of PMS, except for a 
shift of the peak position, which is closer to the tetrahc- 
dral position, as well as for a more pronounced presence 
of loops. Interestingly, when <f> is increased, local tetrahc- 
dral order is lost, as for the reported <f> = 0.56-curve, and 
P(9) tends to the N max distribution, confirming that the 
organization of the bonds becomes increasingly random. 



E. Liquid-Gas unstable region 

Although we do not perform accurate phase coexis- 
tence studies, we can provide a rough estimate of the 
location of the liquid-gas phase separation region and of 
the associated critical point by a combined check of the 
time evolution of the large length-scale structural prop- 
erties, such as S(q — > 0), as well as of the pressure be- 
haviour along different isotherms. The lowest investi- 
gated <p where phase separation can be ruled out at all 
studied temperatures is cf> = 0.25. For lower cf>, we de- 
tect an increase of S(q — » 0) as well as the development 
of maxima and minima in the behaviour of P((f>) (com- 
ing respectively from small and large <fi). We define, for 
each T, the location in of such maxima and minima as 
the spinodal points. We note that at T = 0.1 the pres- 
sure dependence on <fi does not show any loop, while at 
T = 0.095 a small instability region seems to be present 
(within our numerical resolution) close to <f> ~ 0.10. 
These two temperatures thus bracket the critical temper- 
ature T c . Hence, we estimate T c = 0.095 ± 0.005, while 
4> c = 0.10 ± 0.025. A more accurate evaluation of the 
critical parameters, with appropriate techniques should 
be undertaken in future studies. We also note that the 
value of <p at which no critical fluctuations are observed 
in the entire investigated T-range is in close agreement 
both with the -/V max model [TtJ| with coordination number 
four as well as with the PMW[T5] and PMSfHJ models, 
reinforcing the view that the valency is the control pa- 
rameter for the location of the liquid-gas spinodal |11|. 
under single-bond conditions. 



IV. A GLOBAL LOOK AT THE PHASE 
DIAGRAM 

To summarize the static properties and to provide a 
coherent picture of the system behaviour, we discuss here 
the phase diagram of the model, in combination with 
other relevant thermodynamic and connectivity loci. 

Fig. [lO] shows the liquid-gas spinodal, the locus en- 
compassing the region of thermodynamic instability. It is 
worth emphasizing that such instability region is confined 
to densities much smaller than those observed in spher- 
ical attractive potentials. Liquid densities are roughly 
half of the typical value of normal liquids, highlighting 
the large empty spaces which characterize the formation 
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FIG. 10: Static phase diagram in the (</>, T) plane of the stud- 
ied spherical binary model. The reported lines are: gas-liquid 
spinodal (circles), percolation line (x-symbols), locus of CV 
maxima (diamonds), iso-p^ lines at selected values (triangles) 
and energy minima (crosses). Lines are guides to the eye. 
Vertical dashed lines show the liquid boundary of the spin- 
odal line (<f> ~ 0.235) and the boundary (<f) ~ 0.40) between 
the optimal network region, where bonds are primarily direc- 
tional, and that of a disordered network with defects, where 
directionality is progressively lost. 



of a tetrahedral network. At the same time, the small-0 
location of the liquid branch of the spinodal shows that, 
in this model, it is possible to find a region of interme- 
diate densities where, on lowering T, the system remains 
homogeneous down to very small T, a feature which is 
not possible with spherically symmetric attractive poten- 
tials. Moreover, the liquid branch of the spinodal termi- 
nates in a region which is compatible with the window of 
stability of the diamond crvstal[50ll5Tj. This topology of 
the liquid-gas coexistence is typical of models with direc- 
tional interactions, and it is here confirmed through the 
use of our spherical binary mixture. 

To provide an indication of the degree of bonding ob- 
served in different regions of the phase diagram, Fig. [10] 
also shows constant bond-probability lines (iso-p;, lines), 
i.e. lines of constant number of bonds per particles. Here 
a bond is defined as a floating bond with potential en- 
ergy —2uq. The lines have positive slope, signaling that 
on increasing density bond formation becomes favored. 
Only at very large <p an d very large pb (only visible for 
Pb = 0.95 in the figure for the highest density point), a 
sharp change of slope (flat, then negative at even larger 
Pb) is found, associated to the disruptive effect of increas- 
ing density beyond the values for optimal network forma- 
tion. 



FIG. 11: Static phase diagram in the ((f), T) plane of the stud- 
ied spherical binary model. In addition to lines reported in 
Fig. 10 i.e. gas-liquid spinodal (circles), locus of CV maxima 



Fig. 10 also shows the percolation line and the locus 
of Cy maxima. State points on the right-hand side of 
the percolation line are characterized by the presence of 
an infinite cluster in > 50% of the configurations. This 



(diamonds) and of energy minima (crosses), several isobaric 
paths are shown (triangles). Moreover, the locus of Cp max- 
ima (squares) is also reported, merging into the spinodal line 
at the critical point. 



definition of percolation locus is strictly a geometric mea- 
sure and does not provide any information on the lifetime 
of the spanning cluster [5]. The percolation locus coin- 
cides, within numerical resolution, with the pb — 0.45 
locus, suggesting that a constant fraction of bonds is re- 
quested to percolate independently of 4>. The pb value 
at percolation, p v h erc ~ 0.45, is slightly larger than the 
one for random bond percolation on a diamond lattice, 
known to be pb — 0.388[52 , an effect which can be 
attributed to the disordered distribution of particles in 
the fluid. The percolation locus is located above the 
critical point, confirming that the development of long 
range correlated critical fluctuations requires as a pre- 
requisite the existence of a spanning network of bonded 
particles [53]. The locus of C v maxima is located well in- 
side the percolation region and, by comparing with the 
iso-pb curves, it takes place when the bond probability is 
ss 0.75. Comparing with published data of models with 
different valency [39 | [41 ] 142] it appears that the Cy max- 
ima line progressively moves to larger and larger pb values 
with increasing valence. The trend is thus opposite to the 
one characterizing the valence dependence of p^ erc . Fi- 
nally, Fig. 10 also shows the locus of potential energy 
minima along isotherms, i.e. for which (dE / '04>)t = 0, 
a line marking the crossover from an unconstrained net- 
work of bonds to a state in which bonding can not be 
any longer optimized due to volume constraints. 

Additional information on the T and p dependence of 



the properties of the model are reported in Fig. 11 This 
shows, in addition to some of the lines reported in the 
previous figure, the location of selected isobars with the 
pressure P varying by more than two orders of magni- 
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ude. From the constant-P paths, we evaluate the en- 
thalpy H = E + PV, and we extract the behaviour of 
the constant-pressure specific heat Cp = (dH/dT) p. We 
find the presence of clear maxima also in Cp, whose locus 
is also reported in Fig. [IT] The line of Cp maxima is an 



example of a Widom line [54, 55 , which starts, by defini- 
tion, from the critical point. Its calculation therefore also 
confirms the previously discussed estimate of the critical 
point within our numerical resolution. Finally we observe 
that, in the investigated T region, there is no evidence of 
anomalies in the density (constant P density maxima), 
that are commonly found in water models 1561 157] , 



V. CONCLUSIONS 

In this article we have introduced an opportunely de- 
signed, spherical binary mixture model, which is able to 
generate a fluid with a effective coordination number of 
four. Differently from previously studied four-coordinate 
model[inj EJ H2 US E], the particle-particle interac- 
tions are spherical, a feature which is important in order 
to allow a future comparison with theoretical approaches. 

We have shown that the static phase diagram of the 
system is very similar to the one reported for patchy 
models [12l [13j [HI [58], as well as for more sophisticated 
models of network forming liquids [57J [S3]. Indeed, the 
unstable phase-separating region of the system is con- 
fined to low packing fraction (f> and small T, consistently 
with previous studies and Wertheim theory calculations 
[11) for four-coordinated particles. This reinforces the 
statement, that, from the point of view of suppressing 
phase separation, the arrangement of the sticky points 
onto the particle surface is not qualitatively, but only 
quantitatively relevant. 

The unstable gas- liquid region is followed at larger den- 
sities by an optimal network region, i.e. a window of den- 
sities where the system can be equilibrated down to very 
low temperatures observing a progressive formation of a 
tetrahedral network of bonds. In this region, the system 
almost reaches the fully bonded configuration, i.e. the 
disordered ground state of the system and hence a fur- 



ther lowering of the temperature would not produce any 
significant structural change. For even larger densities, 
we observe a destruction of the tetrahedral bonding, in- 
duced by the packing constraints. A signature of this 
effect is observed in both the structure factor of the sys- 
tem, which crosses from the tetrahedral-network form to 
the hard-sphere form as well as in the progressive break- 
ing of the bonds on isothermal increase of the density. 

It is also interesting to discuss the relative location 
of the percolation locus, the locus of specific heat max- 
ima and the liquid-gas spinodal. In the present model, 
the percolation line provides the first indication of the 
clustering process upon lowering T. Successive cooling 
brings to the presence of a line of Cy maxima, and finally 
of the spinodal line. The constant-volume specific heat 
maxima line appears to be a characteristic of all bonded 
systems, since it is observed also in the case of valence 
two [41] where no percolation is present. Its relative loca- 
tion has an opposite dependence on the bond probability 
as compared to the percolation line, when studied as a 
function of the valence. This suggests that, for values 
of the valence larger than the one studied here, the line 
of Cy maxima may become buried in the region where 
equilibration is not feasible any longer, i.e. in the glass 
state. This could be the reason why it is not observed in 
standard spherical models, such as a simple square- well 
potential. 

Finally, we note that the possibility of generating the 
complex pattern characteristic of tetrahedral liquids with 
a simple binary mixture with pair-wise additive square- 
well interactions opens the possibility of applying the 
tools of modern liquid theory (which are nowadays par- 
ticularly accurate for spherical potentials) to the present 
case. It will also be possible, using theoretical and/or nu- 
merical structure factors investigate the dynamics of the 
present system using the formalism of the Mode Coupling 
Theory for the glass transition. Work in this direction is 
underway and a companion paper on the dynamics will 
appear shortly. 
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